// Test my algorithm for dmcs model streaming. We will just simulate
// the belief states with number, and the joinable relationship is
// randomly generated. Then we run our algorithm and compare it with
// the naive approach which enumerate all possible combinations of
// models. Note that when bring this to the real dmcs setting, one has
// to deal with network communication and requerying the neighbors,
// which restarts the whole computation at the neighbor. This is a
// trade off between time and space. (Time is infinite but space is
// bounded).

// to compile:
// g++ -I/opt/local/include /opt/local/lib/libboost_program_options-mt.a dmcs-models-streaming.cpp


#include <algorithm>
#include <fstream>
#include <iostream>
#include <list>
#include <map>
#include <set>
#include <vector>

#include <boost/graph/adjacency_list.hpp>
#include <boost/program_options.hpp>
#include <boost/shared_ptr.hpp>

#define DEBUG
#define MAX_CACHE 20

#define HELP        "help"
#define SYSTEM_SIZE "system-size"
#define LBOUND_ANS  "lower-bound"
#define UBOUND_ANS  "upper-bound"
#define DENSITY     "density"
#define INPUT_GRAPH "graph"

#ifdef DEBUG
#define PREFIX      "prefix"

std::ofstream file_joint_graph;
std::ofstream file_answer_naive;
std::ofstream file_answer_sophisticated;

std::string filename_joint_graph;
std::string filename_answer_naive;
std::string filename_answer_sophisticated;
std::string filename_input_joint_graph;

std::string prefix;
#endif


typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::bidirectionalS> Graph;
typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
typedef boost::graph_traits<Graph>::vertex_iterator   VertexIter;
typedef boost::graph_traits<Graph>::edge_descriptor   Edge;
typedef boost::graph_traits<Graph>::edge_iterator     EdgeIter;

typedef boost::shared_ptr<Graph> GraphPtr;

// a Node simulates an answer (ans_id) from a context (ctx_id)
struct Node
{
  std::size_t ctx_id;
  std::size_t ans_id; // 0 means EOF

  Node(std::size_t c, std::size_t a)
    : ctx_id(c), ans_id(a)
  { }
};

bool
operator== (const Node& n1, const Node& n2)
{
  return ((n1.ctx_id == n2.ctx_id) && (n1.ans_id == n2.ans_id));  
}

typedef boost::shared_ptr<Node> NodePtr;

// a NodeVec simulates a stream of answers from a context
typedef std::vector<NodePtr> NodeVec;
typedef boost::shared_ptr<NodeVec> NodeVecPtr;

typedef std::vector<NodeVec::iterator> NodeVecIter;
typedef boost::shared_ptr<NodeVecIter> NodeVecIterPtr;

// a NodeVec2 simulates the situation when we combine
// models/beliefstates at a local context, taking models from
// neighbors into account. We take one element from each NodeVec into
// account to make it similar to the stream of models.
typedef std::vector<NodeVecPtr> NodeVec2;
typedef boost::shared_ptr<NodeVec2> NodeVec2Ptr;

typedef std::list<NodePtr> NodeList;
typedef boost::shared_ptr<NodeList> NodeListPtr;
typedef std::list<NodeListPtr> NodeList2;
typedef boost::shared_ptr<NodeList2> NodeList2Ptr;

// lazily store the graph as a multimap of (duplicate) pairs of nodes
typedef std::pair<NodePtr, NodePtr> NEdge;
typedef std::multimap<NodePtr, NodePtr> NGraph;
typedef boost::shared_ptr<NGraph> NGraphPtr;

typedef std::list<std::size_t> IntList;
typedef boost::shared_ptr<IntList> IntListPtr;
typedef std::pair<std::size_t, IntListPtr> OneCliqueMask;
typedef std::map<std::size_t, IntListPtr> CliqueMask;
typedef boost::shared_ptr<CliqueMask> CliqueMaskPtr;


CliqueMaskPtr cliquemask(new CliqueMask);

std::size_t system_size;
std::size_t lower_bound;
std::size_t upper_bound;
std::size_t density;


NodeList2Ptr cache(new NodeList2);


NodeVec2Ptr mcs(new NodeVec2);
NGraphPtr joint_graph(new NGraph);

void
find_max_clique(const NGraphPtr joint_graph, NodeListPtr subg, NodeListPtr cand, 
		NodeListPtr clique, NodeList2Ptr maximum_cliques, std::size_t& max_size);


bool
equal(const NodeList& node_list1, const NodeList& node_list2)
{
  assert ( node_list1.size() == node_list2.size() );
  
  NodeList::const_iterator it1 = node_list1.begin();
  NodeList::const_iterator it2 = node_list2.begin();

  for (; it1 != node_list1.end(); ++it1, ++it2)
    {
      if (!(*it1 == *it2))
	{
	  return false;
	}
    }

  return true;
}


bool
cached(const NodeList2Ptr cache, NodeListPtr new_answer)
{
  for (NodeList2::const_iterator it = cache->begin(); it != cache->end(); ++it)
    {
      NodeList old_answer = **it;
      if (equal(old_answer, *new_answer))
	{
	  return true;
	}
    }
  return false;
}


std::ostream&
operator<< (std::ostream& os, const Node& node)
{
  os << "(" << node.ctx_id << ", " << node.ans_id << ")";
  return os;
}


std::ostream&
operator<< (std::ostream& os, const NodeVec& node_vec)
{
  for (NodeVec::const_iterator it = node_vec.begin(); it != node_vec.end(); ++it)
    {
      os << **it << "  ";
    }
  return os;
}


std::ostream&
operator<< (std::ostream& os, const NodeList& node_list)
{
  for (NodeList::const_iterator it = node_list.begin(); it != node_list.end(); ++it)
    {
      //os << "(" << *it << ", " << **it << ")  ";
      os << **it << " ";
    }
  return os;
}

std::ostream&
operator<< (std::ostream& os, const NodeList2& node_list2)
{
  for (NodeList2::const_iterator it = node_list2.begin(); it != node_list2.end(); ++it)
    {
      os << **it << std::endl;
    }
  return os;
}


std::ostream&
operator<< (std::ostream& os, const NodeVec2& node_vec2)
{
  for (NodeVec2::const_iterator it = node_vec2.begin(); it != node_vec2.end(); ++it)
    {
      os << **it << std::endl;
    }
  return os;
}


std::ostream&
operator<< (std::ostream& os, const NGraph& joint_graph)
{
  for (NGraph::const_iterator it = joint_graph.begin(); it != joint_graph.end(); ++it)
    {
      NodePtr v = it->first;
      NodePtr u = it->second;
      os << v->ctx_id << " " << v->ans_id << " " << u->ctx_id << " " << u->ans_id << std::endl;
    }
  return os;
}


std::ostream&
operator<< (std::ostream& os, const IntList& int_list)
{
  for (IntList::const_iterator it = int_list.begin(); it != int_list.end(); ++it)
    {
      os << *it << " ";
    }
  return os;
}


void
generate_relationship(const NodeVec2Ptr mcs, NGraphPtr joint_graph, std::size_t density, std::size_t ci, std::size_t cj)
{
  const NodeVecPtr context_i = (*mcs)[ci];
  const NodeVecPtr context_j = (*mcs)[cj];

  NodeVec::const_iterator end_i = context_i->end();
  --end_i;
  NodeVec::const_iterator end_j = context_j->end();
  --end_j;

  for (NodeVec::const_iterator it_i = context_i->begin(); it_i != end_i; ++it_i)
    {
      for (NodeVec::const_iterator it_j = context_j->begin(); it_j != end_j; ++it_j)
	{
	  // decide whether the 2 answers are joinable with a probability of [density/10]
	  std::size_t flip = rand() % 10;
	  if (flip < density)
	    {
	      // insert to the map
	      NEdge e1(*it_i, *it_j);
	      NEdge e2(*it_j, *it_i);
	      joint_graph->insert(e1);
	      joint_graph->insert(e2);
	    }
	}
    }
}


// Check whether a new answer is joinable with all of the chosen
// answers (from other contexts) which are stored in "result"
bool
joinable(const NodeVecPtr result, const NodePtr it, const NGraphPtr joint_graph)
{
  NGraph::const_iterator lbound = joint_graph->lower_bound(it);
  NGraph::const_iterator ubound = joint_graph->upper_bound(it);

  for (NodeVec::const_iterator jt = result->begin(); jt != result->end(); ++jt)
    {
      bool ok = false;
      for (NGraph::const_iterator gt = lbound; gt != ubound; ++gt)
	{
	  const NEdge e = *gt;
	  if (e.second == *jt)
	    {
	      ok = true;
	      break;
	    }
	}
      if (!ok)
	{
	  return false;
	}
    }

  return true;
}


// just a backtracking, exhaustive search to have the correct number
// of answers for later comparison with our algorithm
void
find_answer_naive(const NodeVec2Ptr mcs, const NGraphPtr joint_graph, std::size_t& no_ans_naive, NodeVecPtr result, std::size_t depth)
{
  // 
  NodeVecPtr context = (*mcs)[depth];

  NodeVec::const_iterator end = context->end();
  --end;
  for (NodeVec::const_iterator it = context->begin(); it != end; ++it)
    {
      if (joinable(result, *it, joint_graph))
	{
	  result->push_back(*it);
	  if (depth+1 == mcs->size())
	    {
	      // got an answer
	      no_ans_naive++;

#ifdef DEBUG
	      file_answer_naive << *result << std::endl;
#endif
	    }
	  else
	    {
	      find_answer_naive(mcs, joint_graph, no_ans_naive, result, depth+1);
	    }
	  result->pop_back();
	}
    }
}


void
prepare_filenames()
{
  filename_joint_graph = prefix + "-jg.txt";
  std::cerr << "filename_joint_graph = " << filename_joint_graph << std::endl;
  filename_answer_naive = prefix + "-ans-n.txt";
  filename_answer_sophisticated = prefix + "-ans-s.txt";

  file_answer_naive.open(filename_answer_naive.c_str());
  file_answer_sophisticated.open(filename_answer_sophisticated.c_str());
}

void
init()
{
  // initialize clique masks
  std::size_t end = 1 << system_size;
  for (std::size_t i = 1; i < end; ++i)
    {
      IntListPtr list(new IntList);

      // naive way to count the number of bits of an integer
      for (std::size_t j = 0; j < system_size; ++j)
	{
	  std::size_t check_bit = 1 << j;
	  if (check_bit & i)
	    {
	      list->push_back(0);
	    }
	}

      (*cliquemask)[i] = list;
    }
}


std::size_t
encoded_index(const NodeListPtr node_list)
{
  std::size_t index = 0;

  //  std::cerr << "node_list = " << *node_list << std::endl;

  for (NodeList::const_iterator it = node_list->begin(); it != node_list->end(); ++it)
    {
      NodePtr node = *it;
      std::size_t ctx_id = node->ctx_id;
      std::size_t shift = 1 << ctx_id;
      index |= shift;
    }

  //  std::cerr << "index = " << index << std::endl;

  return index;
}


bool
already_visited(const NodeListPtr node_list)
{
  std::size_t index = encoded_index(node_list);
  IntListPtr mask = (*cliquemask)[index];

  assert ( mask->size() == node_list->size() );

  std::cerr << "In already visited: " << std::endl
	    << "mask      = " << *mask << std::endl
	    << "node_list = " << *node_list << std::endl;
  
  IntList::const_iterator it = mask->begin();
  NodeList::const_iterator nt = node_list->begin();

  for (; it != mask->end(); ++it, ++nt)
    {
      if (*it < (*nt)->ans_id)
	{
	  return false;
	}
    }

  return true;
}


void
update_clique_mask(const NodeListPtr node_list)
{
  std::size_t index = encoded_index(node_list);
  IntListPtr mask = (*cliquemask)[index];

  assert ( mask->size() == node_list->size() );

  std::cerr << "clique mask[" << index << "] before updating: " << *mask << std::endl;

  IntList::iterator it = mask->begin();
  NodeList::const_iterator nt = node_list->begin();

  for (; it != mask->end(); ++it, ++nt)
    {
      *it = (*nt)->ans_id;
    }

  std::cerr << "clique mask[" << index << "]  after updating: " << *mask << std::endl;
}


void
close_stuff()
{
#ifdef DEBUG
  file_answer_naive.close();
  file_answer_sophisticated.close();
#endif
}


NodeVec::iterator
getFirstModel(NodeVec2Ptr mcs, std::size_t ctx_id)
{
  NodeVec2::const_iterator it = mcs->begin();
  std::advance(it, ctx_id);
  
  NodeVecPtr ctx = *it;
  return ctx->begin();
}


// collect all neighbors of a node into a set
NodeListPtr
getNeighbors(const NodePtr u, const NGraphPtr joint_graph)
{
  NodeListPtr neighbors(new NodeList);
  NGraph::const_iterator lbound = joint_graph->lower_bound(u);
  NGraph::const_iterator ubound = joint_graph->upper_bound(u);

  for (NGraph::const_iterator it = lbound; it != ubound; ++it)
    {
      NEdge e = *it;
      neighbors->push_back(e.second);
    }

  return neighbors;
}


/*void
process(const NodeVec2Ptr mcs, const NGraphPtr joint_graph, 
	NodeVecIterPtr current_models, std::size_t& no_ans_improve,
	NodeListPtr subg, NodeListPtr orig_cand)
{
  std::cerr << "in process()" << std::endl;
  if (subg->size() == mcs->size()) // we have an answer
    {
      no_ans_improve++;
      std::cerr << "This is an answer:" << std::endl
		<< *subg << std::endl;
    }
  else
    { // got a sub, complete graph
      std::cerr << "Got a sub, complete graph:" << std::endl
		<< *subg << std::endl;

      // the rest is the difference between orig_cand and subg
      NodeListPtr the_rest(new NodeList);
      std::back_insert_iterator<NodeList> insert_the_rest(*the_rest);
      orig_cand->sort();
      subg->sort();
      std::set_difference(orig_cand->begin(), orig_cand->end(), 
			      subg->begin(), subg->end(), insert_the_rest);
      
      assert ( the_rest->size() > 0 );
      
      std::cerr << "The rest = " << *the_rest << std::endl;
      
      // now get next answer from each context in "the_rest"
      bool has_eom = false;
      for (NodeList::iterator it = the_rest->begin(); it != the_rest->end(); ++it)
	{
	  std::cerr << "Checking " << **it << std::endl;
	  NodePtr ctx = *it;
	  std::size_t id = ctx->ctx_id;
	  std::cerr << "id = " << id << std::endl;
	  NodeVec::iterator p_ans = (*current_models)[id];
	  std::cerr << "p_ans " << **p_ans << std::endl;
	  p_ans++;
	  NodePtr ans = *p_ans;
	  
	  std::cerr << "Next answer from " << id << " = " << *ans << std::endl;
	  
	  if (ans->ans_id == 0)
	    {
	      std::cerr << "Reached end of models" << std::endl;
	      has_eom = true;
	    }
	  (*current_models)[id] = p_ans;
	}
      
      std::cerr << "finished getting next models for all contexts from \"the rest\"" << std::endl;

      if (has_eom)
	{
	  // there is some ctx out of the clique that reached end of models
	  
	  // pick a ctx in the clique and get next answer from it
	  bool all_end = true;
	  for (NodeList::const_iterator it = subg->begin(); it != subg->end(); ++it)
	    {
	      NodePtr ctx = *it;
	      std::size_t id = ctx->ctx_id;
	      NodeVec::iterator p_ans = (*current_models)[id];
	      p_ans++;
	      NodePtr ans = *p_ans;
	      
	      std::cerr << "Next answer from " << id << " = " << *ans << std::endl;

	      (*current_models)[id] = p_ans;
	      
	      if (ans->ans_id != 0)
		{
		  all_end = false;
		  break;
		}
	    }
	  
	  if (all_end)
	    {
	      // no more hope with this clique
	      std::cerr << "return when all elements in clique reached eof" << std::endl;
	      return;
	    }

	  // otherwise, reset all contexts that reached EOF
	  std::cerr << "Reset current models of EOF guys to first models" << std::endl;
	  NodeVecIter::const_iterator cm_it = current_models->begin();
	  NodeVec2::const_iterator mcs_it = mcs->begin();
	  for (std::size_t i = 0; mcs_it != mcs->end(); ++i, ++mcs_it)
	    {
	      NodeVec::iterator node_iter = (*current_models)[i];
	      NodePtr node = *node_iter;

	      if (node->ans_id == 0)
		{
		  NodeVecPtr ctx = *mcs_it;
		  NodeVec::iterator first_model = ctx->begin();
		  (*current_models)[i] = first_model;
		}
	    }
	}
      
      // reset cand and orig_cand and continue to try
      
      NodeListPtr new_subg(new NodeList);
      NodeListPtr new_cand(new NodeList);
      // NodeListPtr new_fini(new NodeList);
      NodeListPtr clique (new NodeList);
      NodeListPtr new_orig_cand(new NodeList);
      
      std::cerr << "Reset current_models: ";
      for (NodeVecIter::const_iterator it = current_models->begin(); it != current_models->end(); ++it)
	{
	  new_cand->push_back(**it);
	  new_orig_cand->push_back(**it);
	  std::cerr << ***it << " ";
	}
      std::cerr << std::endl;
      
      std::cerr << "Start finding max clique" << std::endl;
      
      find_max_clique(joint_graph, new_subg, new_cand, clique, new_orig_cand);      
    }
  std::cerr << "return from process" << std::endl;
  }*/



void
find_max_clique(const NGraphPtr joint_graph, NodeListPtr subg, NodeListPtr cand, 
		NodeListPtr clique, NodeList2Ptr maximum_cliques, std::size_t& max_size)
{
  //  std::cerr << "subg = " << *subg << std::endl;
  //std::cerr << "cand = " << *cand << std::endl;
  //std::cerr << "clique = " << *clique << std::endl;

  if (subg->empty())
    {
      // Got a maximal clique, process it. For the moment let's just
      // find 1 maximal clique, later we can backtrack and try all
      // different max clique.
      //std::cerr << "GOT A MAXIMAL CLIQUE: " << *clique << std::endl;

      if (clique->size() >= max_size)
	{
	  if (clique->size() > max_size)
	    {
	      maximum_cliques->clear();
	      max_size = clique->size();
	      //std::cerr << "reset to new max_size = " << max_size << std::endl;
	    }
	  
	  NodeListPtr new_clique(new NodeList);
	  for (NodeList::const_iterator it = clique->begin(); it != clique->end(); ++it)
	    {
	      //NodePtr c(new Node(**it));
	      new_clique->push_back(*it);
	    }

	  maximum_cliques->push_back(new_clique);
	}
    }
  else
    {
      // choose a pivot vertex u in cand \cup fini such that |cand \cap N(u)| ~~> max
      std::size_t max = 0;
      NodePtr pivot;
      NodeListPtr intersection;
      NodeListPtr neighbors;
      
      //std::cerr << "In find_max_clique. " << std::endl;
	//		<< "cand = " << *cand << std::endl;
      
      cand->sort();
      
      for (NodeList::const_iterator it = cand->begin(); it != cand->end(); ++it)
	{
	  //std::cerr << "checking for pivot, it = " << **it << ".";
	  NodeListPtr tmp_neighbors = getNeighbors(*it, joint_graph);
	  //std::cerr << "neighbors = " << *tmp_neighbors << std::endl;
	  tmp_neighbors->sort();
	  cand->sort();
	  
	  NodeListPtr tmp_intersection(new NodeList);
	  
	  std::back_insert_iterator<NodeList> insert_tmp_intersection(*tmp_intersection);
	  std::set_intersection(cand->begin(), cand->end(), tmp_neighbors->begin(),
				tmp_neighbors->end(), insert_tmp_intersection);
	  
	  //std::cerr << "intersection with cand = " << *tmp_intersection << std::endl;
	  
	  if (max < tmp_intersection->size() + 1)
	    {
	      max = tmp_intersection->size();
	      pivot = *it;
	      intersection = tmp_intersection;
	      neighbors = tmp_neighbors;
	    }
	}
      
#ifdef DEBUG
      //std::cerr << "Pivot node          = " << *pivot << std::endl
      //	<< "Maximal clique size = " << max << std::endl;
#endif
      
      const NodeListPtr difference(new NodeList);
      if (max > 0)
	{
	  intersection->sort();
	  cand->sort();
	  // now we have the pivot node
	  std::back_insert_iterator<NodeList> insert_difference(*difference);
	  std::set_difference(cand->begin(), cand->end(), intersection->begin(), intersection->end(), insert_difference);
	}
      else // intersection is pointing to NULL now, difference is a copy of cand
	{
	  for (NodeList::const_iterator it = cand->begin(); it != cand->end(); ++it)
	    {
	      difference->push_back(*it);
	    }
	}

      for (NodeList::const_iterator it = difference->begin(); it != difference->end(); ++it)
	{
	  clique->push_front(*it);

	  NodeListPtr new_subg(new NodeList);
	  std::back_insert_iterator<NodeList> insert_new_subg(*new_subg);
	  NodeListPtr new_cand(new NodeList);
	  std::back_insert_iterator<NodeList> insert_new_cand(*new_cand);
	  
	  subg->sort();
	  cand->sort();

	  NodeListPtr chosen_neighbors = getNeighbors(*it, joint_graph);
	  
	  std::set_intersection(subg->begin(), subg->end(), chosen_neighbors->begin(), chosen_neighbors->end(), insert_new_subg);
	  std::set_intersection(cand->begin(), cand->end(), chosen_neighbors->begin(), chosen_neighbors->end(), insert_new_cand);
	  
	  subg->push_front(*it);
	  find_max_clique(joint_graph,new_subg, new_cand, clique, maximum_cliques, max_size);

	  cand->remove(*it);
	  clique->pop_front();
	}
    }
}


void
find_answer_improve(const NodeVec2Ptr mcs, const NGraphPtr joint_graph, 
		    NodeVecIterPtr current_models, std::size_t& no_ans_improve, bool reseted)
{
  NodeListPtr subg(new NodeList);
  NodeListPtr cand(new NodeList);
  NodeListPtr orig_cand(new NodeList);
  NodeListPtr clique(new NodeList);
  NodeList2Ptr maximum_cliques(new NodeList2);
  std::size_t max_size = 0;

  for (NodeVecIter::const_iterator it = current_models->begin(); it != current_models->end(); ++it)
    {
      cand->push_back(**it);
      orig_cand->push_back(**it);
      subg->push_back(**it);
    }
  orig_cand->sort();

  std::cerr << "current graph = " << *cand << std::endl;

  find_max_clique(joint_graph, subg, cand, clique, maximum_cliques, max_size);

  std::cerr << "MAXIMUM CLIQUES: " << std::endl
	    << *maximum_cliques << std::endl;

  //  std::cerr << "PRESS ENTER TO CONTINUE..." << std::endl;
  //char dump_c;
  //std::cin >> dump_c;

  std::cerr << "max_size = " << max_size << std::endl;
  std::cerr << "mcs->size = " << mcs->size() << std::endl;
  if (max_size == mcs->size())
    {
      // all are answers
      std::cerr << "FOUND SOME ANSWER. PRESS ENTER TO CONTINUE..." << std::endl;
      //      no_ans_improve += maximum_cliques->size();


      for (NodeList2::const_iterator it = maximum_cliques->begin(); it != maximum_cliques->end(); ++it)
	{
	  if (!cached(cache, *it))
	    {
	      cache->push_back(*it);
	      no_ans_improve++;
	    }

#ifdef DEBUG
	  file_answer_sophisticated << **it << std::endl;
	}
#endif

      // find a node to flip
      bool all_end = true;

      bool has_eom = false;
      NodeVecIterPtr new_current_models(new NodeVecIter);
      for (NodeVecIter::const_iterator it = current_models->begin(); it != current_models->end(); ++it)
	{
	  new_current_models->push_back(*it);
	}

      for (std::size_t i = 0; i < system_size; ++i)
	{
	  NodeVec::iterator p_ans = (*current_models)[i];
	  NodeVec::iterator next_ans = p_ans;
	  next_ans++;
	  NodePtr ans = *next_ans;

	  if (ans->ans_id != 0)
	    {
	      all_end = false;
	      (*new_current_models)[i] = next_ans;

	      std::cerr << "go to explore new answers after finding one" << std::endl;
	      find_answer_improve(mcs, joint_graph, new_current_models, no_ans_improve, has_eom);

	      (*new_current_models)[i] = p_ans;
	    }
	}

      if (all_end)
	{
	  std::cerr << "return when all elements in clique reached eof" << std::endl;
	  return;
	}
    }
  else
    {
      // for each maximum clique (sub-complete graph), try to flip the outside nodes
      for (NodeList2::iterator it = maximum_cliques->begin(); it != maximum_cliques->end(); ++it)
	{

	  NodeVecIterPtr new_current_models(new NodeVecIter);
	  for (NodeVecIter::const_iterator jt = current_models->begin(); jt != current_models->end(); ++jt)
	    {
	      new_current_models->push_back(*jt);
	    }

	  NodeListPtr nlp = *it;

	  std::cerr << "pick        " << *nlp << std::endl;
	  std::cerr << "orig_cand = " << *orig_cand << std::endl;

	  if (reseted)
	    {
	      if (already_visited(nlp))
		{
		  std::cerr << "was reseted before and it was already visited! Now RETURN" << std::endl;
		  return;
		}
	      
	      std::cerr << "update clique mask" << std::endl;
	      update_clique_mask(nlp);
	    }

	  NodeListPtr the_rest(new NodeList);
	  std::back_insert_iterator<NodeList> insert_the_rest(*the_rest);

	  nlp->sort();
	  std::set_difference(orig_cand->begin(), orig_cand->end(), nlp->begin(), nlp->end(), insert_the_rest);
	  
	  assert ( the_rest->size() > 0 );

	  // now get next answer from each context in "the_rest"
	  std::cerr << "the_rest = " << *the_rest << std::endl;
	  bool has_eom = false;
	  for (NodeList::iterator jt = the_rest->begin(); jt != the_rest->end(); ++jt)
	    {
	      //std::cerr << "Checking " << **jt << std::endl;
	      NodePtr ctx = *jt;
	      std::size_t id = ctx->ctx_id;
	      //std::cerr << "id = " << id << std::endl;
	      NodeVec::iterator p_ans = (*current_models)[id];
	      //std::cerr << "p_ans " << **p_ans << std::endl;
	      p_ans++;
	      NodePtr ans = *p_ans;
	      
	      //std::cerr << "Next answer from " << id << " = " << *ans << std::endl;
	      
	      if (ans->ans_id == 0)
		{
		  std::cerr << "Reached end of models" << std::endl;
		  // new strategy: if a node outside of the clique
		  // reaches EOF, then all outside nodes will be reset
		  // later, hence we can break here.
		  has_eom = true;
		  break;
		}
	      (*new_current_models)[id] = p_ans;
	    }
      
	  //std::cerr << "finished getting next models for all contexts from \"the rest\"" << std::endl;

	  if (has_eom)
	    {
	      // there is some ctx out of the clique that reached end of models
	      std::cerr << "saw EOF in the_rest" << std::endl;

	      // reset all outside nodes
	      std::cerr << "reset outside nodes" << std::endl;
	      for (NodeList::iterator jt = the_rest->begin(); jt != the_rest->end(); ++jt)
		{
		  //std::cerr << "Checking " << **jt << std::endl;
		  NodePtr ctx = *jt;
		  std::size_t id = ctx->ctx_id;
		  //std::cerr << "id = " << id << std::endl;

		  NodeVecPtr models = (*mcs)[id];
		  NodeVec::iterator first_model = models->begin();
		  (*new_current_models)[id] = first_model;
		}
	      
	      // pick a ctx in the clique, get next answer from it,
	      // and go into the recursive call.
	      bool all_end = true;
	      for (NodeList::const_iterator it = nlp->begin(); it != nlp->end(); ++it)
		{
		  NodePtr ctx = *it;
		  std::cerr << "ctx = " << *ctx << std::endl;
		  std::size_t id = ctx->ctx_id;
		  NodeVec::iterator p_ans = (*current_models)[id];
		  NodeVec::iterator next_ans = p_ans;
		  std::cerr << "id = " << id << ", p_ans" << **p_ans << std::endl;
		  next_ans++;
		  NodePtr ans = *next_ans;
		  
		  //std::cerr << "Next answer from " << id << " = " << *ans << std::endl;
		  

		  
		  if (ans->ans_id != 0)
		    {
		      (*new_current_models)[id] = next_ans;

		      std::cerr << "go into the recursive call" << std::endl;
		      find_answer_improve(mcs, joint_graph, new_current_models, no_ans_improve, has_eom);

		      // restore the value of this node
		      (*new_current_models)[id] = p_ans;

		      all_end = false;
		    }
		}
	  
	      if (all_end)
		{
		  // no more hope with this clique
		  std::cerr << "return when all elements in clique reached eof" << std::endl;
		  return;
		}
	    }
	  else
	    {
	      // go into the recursive call
	      //	  std::cerr << "picked " << **it << std::endl;
	      std::cerr << "go into the recursive call" << std::endl;
	      find_answer_improve(mcs, joint_graph, new_current_models, no_ans_improve, has_eom);
	    }
	}
      
    }
}


void
generate_input()
{
  srand( time(NULL) );
  std::size_t additional_ans_bound = upper_bound - lower_bound;

  for (std::size_t i = 0; i < system_size; ++i)
    {
      // randomly generate the number of answer between lower_bound
      // and upper_bound
      std::size_t no_ans = 0;
      if (additional_ans_bound > 0)
	{
	  std::size_t additional_ans = rand() % additional_ans_bound;
	  no_ans = lower_bound + additional_ans;
	}
      else
	{
	  no_ans = lower_bound;
	}

      NodeVecPtr context(new NodeVec);
      mcs->push_back(context);

      for (std::size_t j = 1; j < no_ans; ++j)
	{
	  NodePtr node(new Node(i, j));
	  context->push_back(node);
	}
      NodePtr node(new Node(i, 0)); // simulation for EOF
      context->push_back(node);
    }

  std::cerr << "The MCS:" << std::endl
	    << *mcs << std::endl;

  // Now generate the "joint relationshp" graph
  for (std::size_t ci = 0; ci < system_size - 1; ++ci)
    {
      for (std::size_t cj = ci+1; cj < system_size; ++cj)
	{
	  generate_relationship(mcs, joint_graph, density, ci, cj);
	}
    }

#ifdef DEBUG
  std::cerr << "Write joint graph to " << filename_joint_graph << std::endl;
  file_joint_graph.open(filename_joint_graph.c_str());
  file_joint_graph << mcs->size() << std::endl;
  for (NodeVec2::const_iterator it = mcs->begin(); it != mcs->end(); ++it)
    {
      NodeVecPtr node = *it;
      file_joint_graph << node->size()-1 << std::endl;
    }
  file_joint_graph << *joint_graph << std::endl;
  file_joint_graph.close();
#endif
}


void
manual_input()
{
  NodeVecPtr context0(new NodeVec);
  NodeVecPtr context1(new NodeVec);
  NodeVecPtr context2(new NodeVec);
  mcs->push_back(context0);
  mcs->push_back(context1);
  mcs->push_back(context2);

  system_size = mcs->size();


  NodePtr node00(new Node(0, 0));
  NodePtr node01(new Node(0, 1));
  NodePtr node02(new Node(0, 2));
  NodePtr node03(new Node(0, 3));
  context0->push_back(node01);
  context0->push_back(node02);
  context0->push_back(node03);
  context0->push_back(node00);

  NodePtr node10(new Node(1, 0));
  NodePtr node11(new Node(1, 1));
  NodePtr node12(new Node(1, 2));
  NodePtr node13(new Node(1, 3));
  context1->push_back(node11);
  context1->push_back(node12);
  context1->push_back(node13);
  context1->push_back(node10);
  
  NodePtr node20(new Node(2, 0));
  NodePtr node21(new Node(2, 1));
  NodePtr node22(new Node(2, 2));
  context2->push_back(node21);
  context2->push_back(node22);
  context2->push_back(node20);


  NEdge e1(node01, node11);
  NEdge e2(node01, node12);
  NEdge e3(node01, node13);
  NEdge e4(node01, node21);
  NEdge e5(node01, node22);
  NEdge e6(node02, node12);
  NEdge e7(node02, node13);
  NEdge e8(node02, node21);
  NEdge e9(node02, node22);
  NEdge e10(node03, node11);
  NEdge e11(node03, node12);
  NEdge e12(node03, node13);
  NEdge e13(node03, node21);
  NEdge e14(node03, node22);
  NEdge e15(node21, node01);
  NEdge e16(node21, node02);
  NEdge e17(node21, node03);
  NEdge e18(node21, node12);
  NEdge e19(node21, node13);
  NEdge e20(node11, node01);
  NEdge e21(node11, node03);
  NEdge e22(node12, node01);
  NEdge e23(node12, node02);
  NEdge e24(node12, node03);
  NEdge e25(node12, node21);
  NEdge e26(node12, node22);
  NEdge e27(node13, node01);
  NEdge e28(node13, node02);
  NEdge e29(node13, node03);
  NEdge e30(node13, node21);
  NEdge e31(node13, node22);
  NEdge e32(node22, node01);
  NEdge e33(node22, node02);
  NEdge e34(node22, node03);
  NEdge e35(node22, node12);
  NEdge e36(node22, node13);

/*  NEdge e0(node11, node21);
  joint_graph->insert(e0);

*/

/*  NEdge e0(node01, node11);
  NEdge e1(node01, node12);
  NEdge e2(node01, node13);
  NEdge e3(node01, node21);
  NEdge e4(node01, node22);
  NEdge e5(node02, node12);
  NEdge e6(node02, node13);
  NEdge e7(node02, node21);
  NEdge e8(node02, node22);
  /*NEdge e9(node03, node11);
  NEdge e10(node03, node12);
  NEdge e11(node03, node13);
  NEdge e12(node03, node21);
  NEdge e13(node03, node22);
  NEdge e14(node21, node01);
  NEdge e15(node21, node02);
  //NEdge e16(node21, node03);
  NEdge e17(node21, node12);
  NEdge e18(node21, node13);
  NEdge e19(node11, node01);
  //NEdge e20(node11, node03);
  NEdge e21(node12, node01);
  NEdge e22(node12, node02);
  //NEdge e23(node12, node03);
  NEdge e24(node12, node21);
  NEdge e25(node12, node22);
  NEdge e26(node13, node01);
  NEdge e27(node13, node02);
  //NEdge e28(node13, node03);
  NEdge e29(node13, node21);
  NEdge e30(node13, node22);
  NEdge e31(node22, node01);
  NEdge e32(node22, node02);
  //NEdge e33(node22, node03);
  NEdge e34(node22, node12);
  NEdge e35(node22, node13);*/

/*  NEdge e0(node01, node11);
  NEdge e1(node01, node13);
  NEdge e2(node01, node22);
  NEdge e3(node02, node11);
  NEdge e4(node02, node12);
  NEdge e5(node02, node13);
  NEdge e6(node02, node22);
  NEdge e7(node21, node11);
  NEdge e8(node21, node12);
  NEdge e9(node11, node01);
  NEdge e10(node11, node02);
  NEdge e11(node11, node21);
  NEdge e12(node12, node02);
  NEdge e13(node12, node21);
  NEdge e14(node12, node22);
  NEdge e15(node13, node01);
  NEdge e16(node13, node02);
  NEdge e17(node13, node22);
  NEdge e18(node22, node01);
  NEdge e19(node22, node02);
  NEdge e20(node22, node12);
  NEdge e21(node22, node13);*/

  // joint_graph->insert(e0);
  joint_graph->insert(e1);
  joint_graph->insert(e2);
  joint_graph->insert(e3);
  joint_graph->insert(e4);
  joint_graph->insert(e5);
  joint_graph->insert(e6);
  joint_graph->insert(e7);
  joint_graph->insert(e8);
  joint_graph->insert(e9);
  joint_graph->insert(e10);
  joint_graph->insert(e11);
  joint_graph->insert(e12);
  joint_graph->insert(e13);
  joint_graph->insert(e14);
  joint_graph->insert(e15);
  joint_graph->insert(e16);
  joint_graph->insert(e17);
  joint_graph->insert(e18);
  joint_graph->insert(e19);
  joint_graph->insert(e20);
  joint_graph->insert(e21);
  joint_graph->insert(e22);
  joint_graph->insert(e23);
  joint_graph->insert(e24);
  joint_graph->insert(e25);
  joint_graph->insert(e26);
  joint_graph->insert(e27);
  joint_graph->insert(e28);
  joint_graph->insert(e29);
  joint_graph->insert(e30);
  joint_graph->insert(e31);
  joint_graph->insert(e32);
  joint_graph->insert(e33);
  joint_graph->insert(e34);
  joint_graph->insert(e35);
  joint_graph->insert(e36);

  std::cerr << "joint graph: " << *joint_graph << std::endl;
}


NodePtr
node_position(const NodeVec2Ptr mcs, std::size_t ctx_id, std::size_t ans_id)
{
  NodeVec2::const_iterator it = mcs->begin();
  std::advance(it, ctx_id);

  assert (it != mcs->end() );

  NodeVecPtr current_context = *it;
  NodeVec::const_iterator jt = current_context->begin();
  std::advance(jt, ans_id - 1);

  NodePtr node = *jt;

  assert (node->ans_id != 0);

  return node;
}


void
read_joint_graph(const std::string& filename_input_joint_graph)
{
  std::ifstream file_input_joint_graph(filename_input_joint_graph.c_str());
  if (file_input_joint_graph.is_open())
    {
      std::cerr << "READING FROM FILE" << std::endl;
      file_input_joint_graph >> system_size;
      std::cerr << system_size << std::endl;

      for (std::size_t i = 0; i < system_size; ++i)
	{
	  NodeVecPtr context(new NodeVec);
	  mcs->push_back(context);
	  std::size_t answer_size;
	  file_input_joint_graph >> answer_size;
	  std::cerr << answer_size << std::endl;

	  for (std::size_t j = 1; j <= answer_size; ++j)
	    {
	      NodePtr node(new Node(i, j));
	      context->push_back(node);
	    }
	  NodePtr node_null(new Node(i, 0));
	  context->push_back(node_null);
	}

      while (!file_input_joint_graph.eof())
	{
	  std::size_t node1_ctx;
	  std::size_t node1_ans;
	  std::size_t node2_ctx;
	  std::size_t node2_ans;

	  file_input_joint_graph >> node1_ctx;
	  file_input_joint_graph >> node1_ans;
	  file_input_joint_graph >> node2_ctx;
	  file_input_joint_graph >> node2_ans;
	  
	  NodePtr node1 = node_position(mcs, node1_ctx, node1_ans);
	  NodePtr node2 = node_position(mcs, node2_ctx, node2_ans);
	  NEdge e(node1, node2);
	  joint_graph->insert(e);
	}

      std::cerr << *joint_graph << std::endl;
    }
  else
    {
      std::cerr << "File " << filename_input_joint_graph << " not found!" << std::endl;
      exit(1);
    }
}

int main(int argc, char* argv[])
{
  boost::program_options::options_description desc("Allowed options");
  desc.add_options()
    (HELP, "Help and usage message")
    (SYSTEM_SIZE, boost::program_options::value<std::size_t>(&system_size)->default_value(3), "Set number of contexts in the system")
    (LBOUND_ANS, boost::program_options::value<std::size_t>(&lower_bound)->default_value(3), "Set lower bound for the number of answers from each context")
    (UBOUND_ANS, boost::program_options::value<std::size_t>(&upper_bound)->default_value(5), "Set upper bound for the number of answers from each context")
    (DENSITY, boost::program_options::value<std::size_t>(&density)->default_value(8), "Set density of the randomized \"joint relationship\" graph")
    (INPUT_GRAPH, boost::program_options::value<std::string>(&filename_input_joint_graph)->default_value(""), "Set file name to read the joint graph")
#ifdef DEBUG
    (PREFIX, boost::program_options::value<std::string>(&prefix)->default_value("test"), "Set prefix for output file names")
#endif
    ;

  boost::program_options::variables_map vm;
  boost::program_options::store(boost::program_options::parse_command_line(argc, argv, desc), vm);
  boost::program_options::notify(vm);

  if (vm.count(HELP))
    {
      std::cerr << desc << std::endl;
      return 1;
    }

  // Empty systems are not allowed
  assert ( system_size > 0 );

  // We want some realistic bound
  assert ( lower_bound <= upper_bound );

#ifdef DEBUG
  prepare_filenames();
#endif

  if (filename_input_joint_graph.compare("") == 0)
    {
      generate_input();
      //manual_input();
    }
  else
    {
      read_joint_graph(filename_input_joint_graph);
    }
  
  init();

  // Now do the naive enumeration approach just to make sure that we
  // get the correct answers from the graph. This is a recursive
  // procedue.

  std::size_t no_ans_naive = 0;
  NodeVecPtr result(new NodeVec);
  find_answer_naive(mcs, joint_graph, no_ans_naive, result, 0);

  // try models streaming approach. Firstly simulate that we get first
  // model from each context
  NodeVecIterPtr current_models(new NodeVecIter);

  for (NodeVec2::const_iterator it = mcs->begin(); it != mcs->end(); ++it)
    {
      NodeVecPtr ctx = *it;
      NodeVec::iterator first_model = ctx->begin();
      current_models->push_back(first_model);
    }

  std::size_t no_ans_improve = 0;
  find_answer_improve(mcs, joint_graph, current_models, no_ans_improve, false);

  std::cerr << "no_ans_naive = " << no_ans_naive << std::endl;
  std::cerr << "no_ans_improve = " << no_ans_improve << std::endl;
  
  close_stuff();
}
